{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform a 4-terminal calculation with 2 crossed Carbon chains.  \n",
    "Running a two-terminal calculation with TranSiesta is a breeze compared to running $N>2$-electrode calculations. When performing $N>2$-electrode calculations, an endless combination of different applied bias settings become apparent.  \n",
    "This will be reflected in an even more verbose input for TranSiesta to describe all the 4 electrodes, contours, chemical potentials etc.\n",
    "\n",
    "This example will primarily create the geometries, and then you should perform data analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain = sisl.Geometry([[0, 0, 0]], atoms=sisl.Atom[6], lattice=[1.4, 1.4, 11])\n",
    "elec_x = chain.tile(4, axis=0).add_vacuum(11 - 1.4, 1)\n",
    "elec_x.write(\"ELEC_X.fdf\")\n",
    "elec_y = chain.tile(4, axis=1).add_vacuum(11 - 1.4, 0)\n",
    "elec_y.write(\"ELEC_Y.fdf\")\n",
    "chain_x = elec_x.tile(4, axis=0)\n",
    "chain_y = elec_y.tile(4, axis=1)\n",
    "chain_x = chain_x.translate(-chain_x.center(what=\"xyz\"))\n",
    "chain_y = chain_y.translate(-chain_y.center(what=\"xyz\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = chain_x.append(chain_y.translate([0, 0, -chain.cell[2, 2] + 2.1]), 2)\n",
    "# Correct the y-direction vacuum\n",
    "device = device.add_vacuum(chain_y.cell[1, 1] - chain_x.cell[1, 1], 1)\n",
    "device = device.translate(device.center(what=\"cell\"))\n",
    "device.write(\"DEVICE.fdf\")\n",
    "device.write(\"DEVICE.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Run the two electrodes (`RUN_ELEC_X/Y.fdf`).\n",
    "- Run the TranSiesta analyzation step (see fdf flag: `TS.Analyze`) and determine the optimal pivoting scheme used.  \n",
    "  If you are interested you may try to use the worst pivoting scheme and see if it affects the execution time (however this system is very small, so the time difference may be very small).\n",
    "- After analyzation and adding the resulting pivoting scheme to `RUN.fdf`; run the device (`RUN.fdf`).\n",
    "- Try and extract similar data as done in [TB 6](../TB_06/run.ipynb). At least plot one of the DOS quantities.\n",
    "- Extend your DOS plot to be *orbitally resolved* by extracting only subsets of DOS, in this regard also play with the `norm` keyword, try and plot the DOS per $s$, sum of $p$, etc. for the orbitals on the Carbon atoms.\n",
    "\n",
    "  - A file named `siesta.ORB_INDX` has been created by Siesta which contains the orbital information per atom, this should give you access to the indices for extraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbt = sisl.get_sile(\"siesta.TBT.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plot bond (vector) transmissions, below is a skeleton code to do this, look in the `sisl` manual for extraction of *vector transmission*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy = tbt.geometry.xyz[:, :]\n",
    "J1 = # fill in the corresponding code here ()\n",
    "plt.quiver(xy[:, 0], xy[:, 1], J1[:, 0], J1[:, 1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Performing NEGF calculations for $N$-electrodes with an applied bias is *extremely* difficult, and one should first take on this task once the traditional TranSiesta 2-electrode setup is a breeze.\n",
    "\n",
    "One of the main difficulties in performing *good* $N$-electrode calculations is the Poisson solution and the boundary conditions. These are more difficult to impose when having more than 2 electrodes or when having only 1 electrode (only 2 electrode setups are *easy*).\n",
    "\n",
    "- In a 2 terminal device it is obvious how the applied bias is located (a ramp between the two electrodes), however,\n",
    "  when dealing with more than 2 electrodes all electrodes may have a different chemical potential and thus the variations in how to apply the bias becomes infinite.  \n",
    "  Your first task is to read through `RUN.fdf` and figure out which electrode has which chemical potential and draw a small schematic of it.\n",
    "- Calculate the NEGF with a bias of 0.5 eV (please use this command line, for details of options refer to the Siesta manual):\n",
    "\n",
    "      siesta -L no_guess_0.5 -V 0.5:eV RUN.fdf > no_guess_0.5.out\n",
    "      \n",
    "  which applies a bias of 0.5 eV. Read through the output and find the warning which justifies the name `no_guess`.\n",
    "- There should now be 2 files that lets you plot the bias potential profile, `siesta.VH` and `no_guess_0.5.VH`.\n",
    "   Use the below method to plot the plane right between the two carbon chains (*HINT*: calculating the $z$-value at the **centre** between the two chains may be done using a method for the `Geometry` object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_grid(grid, plane_dist=1):\n",
    "    \"\"\" Plot a cut through the grid \"\"\"\n",
    "    z_index = grid.index(plane_dist, 2)or\n",
    "    x, y = np.mgrid[:grid.shape[0], :grid.shape[1]]\n",
    "    dcell = grid.dcell\n",
    "    x, y = x * dcell[0, 0] + y * dcell[1, 0], x * dcell[0, 1] + y * dcell[1, 1]\n",
    "    fig, ax = plt.subplots(1, 1)\n",
    "    im = ax.contourf(x, y, grid.grid[:, :, z_index])\n",
    "    ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]')\n",
    "    ax.set_title('Potential difference [eV]')\n",
    "    # Also plot the atomic coordinates\n",
    "    xyz = grid.geometry.xyz\n",
    "    ax.scatter(xyz[:, 0], xyz[:, 1], 50, 'k', alpha=.6)\n",
    "    fig.colorbar(im);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the two different grids:\n",
    "grid0 = sisl.get_sile(\"siesta.VH\").read_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_guess = sisl.get_sile(\"no_guess_0.5.VH\").read_grid()\n",
    "# Specify the geometry so we can add the atoms to the plot\n",
    "no_guess.set_geometry(device)\n",
    "plot_grid(no_guess - grid0, 1.0)  # replace with the correct z-distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- TranSiesta's method of setting boundary conditions for $N$-electrodes is extremely crude, since it only fixes the potential on the electrodes. Instead of letting TranSiesta apply the boundary conditions, we can provide an external solution to the Poisson problem with *proper* boundary conditions of the electrodes.  \n",
    "  Below is a method to solve the Poisson problem in Python using `pyamg`. It takes quite some time, so be patient.  \n",
    "  You don't need to understand the below script (but I won't hold you back if you want to carefully read it through ;))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the boundary conditions in the unit-cell\n",
    "bc = [[\"dirichlet\", \"dirichlet\"], [\"dirichlet\", \"dirichlet\"], [\"neumann\", \"neumann\"]]\n",
    "\n",
    "# Import the required machinery for solving the boundary conditions\n",
    "# There is also a command-line utility to do this from a siesta.TBT.nc file with\n",
    "# some easier to use command-line options. It isn't fully automated but almost.\n",
    "from sisl_toolbox.transiesta.poisson import solve_poisson\n",
    "\n",
    "device_name = device.copy()\n",
    "\n",
    "# define the electrodes in the device together with their potential\n",
    "elecs = {}\n",
    "# X-left\n",
    "device_name[\"elec-x-1\"] = np.arange(elec_x.na)\n",
    "elecs[\"elec-x-1\"] = 1.0\n",
    "# X-right\n",
    "device_name[\"elec-x-2\"] = np.arange(chain_x.na - elec_x.na, chain_x.na)\n",
    "elecs[\"elec-x-2\"] = 1.0\n",
    "# Y-left\n",
    "device_name[\"elec-y-1\"] = np.arange(chain_x.na, chain_x.na + elec_y.na)\n",
    "elecs[\"elec-y-1\"] = -1.0\n",
    "# Y-right\n",
    "device_name[\"elec-y-2\"] = np.arange(device.na - elec_y.na, device.na)\n",
    "elecs[\"elec-y-2\"] = -1.0\n",
    "\n",
    "# Now solve\n",
    "print(\"Starting solution... Please hold... This takes time! :)\")\n",
    "grid = solve_poisson(\n",
    "    device_name,\n",
    "    grid0.shape,\n",
    "    dtype=np.float32,\n",
    "    tolerance=2e-6,\n",
    "    boundary=bc,\n",
    "    radius=1.7,\n",
    "    **elecs\n",
    ")\n",
    "print(\"Done... storing boundary conditions to disk\")\n",
    "grid.write(\"V.TSV.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- After having solved the Poisson problem, an additional file is created `V.TSV.nc`. It contains the Poisson solution with the appropriate boundary conditions. Now run TranSiesta like this:\n",
    "\n",
    "      siesta -L guess_0.5 -V 0.5:eV -fdf TS.Poisson:V.TSV.nc RUN.fdf > guess_0.5.out\n",
    "      \n",
    "- Once completed, there will be two solutions of the same applied bias `no_guess_0.5` and `guess_0.5`. It is instructive to compare the two:\n",
    "  1. Compare the output of the two runs, is there a difference in the SCF?\n",
    "  2. Plot the potential cut for the `guess_0.5` file as done above for `no_guess_0.5`, how do they differ? Try and plot a few other plane-cuts.\n",
    "  3. Calculate the transmission for both runs using these commands:\n",
    "\n",
    "         tbtrans -L guess_0.5 -V 0.5:eV -D guess_0.5 RUN.fdf > TBT_guess_0.5.out\n",
    "         tbtrans -L no_guess_0.5 -V 0.5:eV -D no_guess_0.5 RUN.fdf > TBT_no_guess_0.5.out\n",
    "      \n",
    "     Plot the two transmissions and density of states, are they the same?\n",
    "- If you have VMD/XCrySDen installed try out the command-line utility `sgrid` to create a cube file containing the potential profile, [this page](https://zerothi.github.io/sisl/scripts/sgrid.html) may be useful. For potential profiles you'll likely prefer to visualize ising `VolumeSlice` for VMD, since that gives how the potential crosse a plane."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
